A Density Matrix Renormalization Group Study of 
Ultrasmall Superconducting Grains 



ON 

on 

On 



(N 



J. Dukelsky 1 and G. Sierra 2 
1 Institute de Estructura de la Materia, C. S. I. C, Madrid, Spain. 
2 Institute/ de Matemdticas y Fisica Fundamental, C.S.I.C., Madrid, Spain. 



We apply the DMRG method to the BCS pairing Hamiltonian which describes ultrasmall super- 
conducting grains. Our version of the DMRG uses the particle (hole) states around the Fermi level 
as the system block (environment). We observe a smooth logarithmic-like crossover between the 
few electron regime and the BCS-bulk regime. 
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The Density Matrix Renormalization Group (DMRG) 
|jj has been applied with great success to a large vari- 
ety of systems in Condensed Matter and Statistical Me- 
chanics (sec Q for a review). Most of these applications 
use the real space formulation although the DMRG can 
also be formulated in momentum space j^]. The start- 
ing point of the DMRG method is the breaking of the 
system into two pieces, the block and the environment, 
which are described by a finite number of states out of 
which one can reconstruct the ground state and the ex- 
cited states of the whole system. A correct choice of the 
block and the environment is therefore crucial for the 
DMRG to work. For example open chains are divided 
into left and right handed pieces linked by a couple of 
sites. The local interaction between the left and right 
pieces of the chains explains the adequacy of the DMRG 
in Id problems. However for more complicated systems, 
or in dimensions higher than one, there are no general 
rules dictating the correct DMRG breaking except the 
nature of the physical problem under study. 

In this letter we shall present a DMRG analysis of the 
ground state (GS) of the BCS Hamiltonian which is used 
to describe the superconducting properties of ultrasmall 
Al grains, discovered by Ralph, Black and Tinkham [Q. 
We shall show that the DMRG gives an accurate approx- 
imation to the exact GS if the block is taken to be the set 
of particles while the environment is taken to be the set 
of holes. This choice does not satisfy the local interaction 
rule, which is so effective in real space, for particles and 
holes are all coupled by the BCS Hamiltonian. Neverthe- 
less the projection of the GS into the particle and hole 
subspaces via density matrices is strongly peaked on a 
small number of states, which explains the applicability 
of the DMRG to this problem. 

The BCS pairing Hamiltonian used for small metallic 
grains is given by J| [TlJ 
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where i,j = 1, 2, . . . , Q label single particle energy levels 
whose energies are given for simplicity by ej = jd, and 
d is the average level spacing which is inversely propor- 
tional to the size of the grain. Cj j0 - are electron destruc- 
tion operators of time reserved states a = ±. Finally fi 
is the chemical potential and A is the BCS coupling con- 
stant, whose appropriate value for the Al grains is 0.224 
H . Given N electrons they can form P Cooper pairs and 
b unpaired states such that N = 2P + b. The Hamilto- 
nian (|l|) decouples the unpaired electrons and hence b is 
a conserved quantity. We shall investigate in this letter 
the GS of even grains (b = 0) and odd grains (b = 1) at 
half filling ( N = f2) using a new version of the DMRG 
and compare our results with the exact Lanczos diago- 
nalization for N = 24 and projected BCS (PBCS) for 
larger N. 

The Hamiltonian ([l]) has two regimes depending on the 
ratio d/A = 2 sinh(l/A)/iV, between the level spacing d 
and the bulk superconducting gap A [|]-|llj . In the weak 
coupling region (d/A » 1), which corresponds to small 
grains or small coupling constants, the system is in a 
regime with strong pairing fluctuations above the Fermi 
sea which lead to logarithmic renormalizations Q . In the 
strong coupling regime {d/A << 1), which corresponds 
to large grains or strong coupling constants, the bulk- 
BCS wave function describes correctly the GS properties. 
In mean field theory the crossover between the weak and 
strong coupling regimes occurs at a critical value of the 
level spacing df which is parity dependent. For even 
grains one has d$ / A = 3.56 while for odd grains c?f / A = 
0.89 IB. 



It is illustrative to consider the case where e 7 - 
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so that the exact GS for b - 
the PBCS state 

\N = ft, b = 0, ej = 0) = 



and N = £1 is given by 
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where |0) is the vacuum of the electron operators and 

Z m ,n = (Ml) 2 C N . M (Cn.m = ^j) is the norm 

of the state (pt) M |0). Notice that (g) is the projection 
of the BCS state exp(P^)|0) into the subspace with n 
electrons. On the other hand for A = and ej = jd the 
GS at half filling is a Fermi sea with all the levels between 
i = l and i = N/2 occupied. This state can be written 
as 



\FS) = U 
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Comparing (g) and (g) it is clear that the PBCS 
state (||) can be derived from the state (||) acting with 
pairs of particle-hole (p-h) creation operators. With 
this aim we split the operator P> as + B\ where 



^ — Ei=n/2+l H,+ H 
HI, 



cj _ . After some algebra one finds 



\n = ft, b = o, ej = o) = Efco ^ I*)* ® 10ft 
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Performing the p-h transformation B —* B\ we deduce 
that (g) can be written as a sum over the tensor product 
of particle and holes states with amplitude ip£, where £ 
is the associated occupation number. Tracing over the 
hole states one obtains from (|J) a density matrix whose 
eigenstates are the particle states \t) v with eigenvalues 
u>e = i>£ [t = 0, 1, . . . , ft/2). Tracing over the particle 
states yield identical results for the hole states. In both 
cases the eigenvalues of the density matrix follow the hy- 
pergeometric distribution we = Cq/ 2 ^/Cn,u/2 which for 
large values of ft becomes a normal distribution centered 
at £1/4 with quadratic deviation a = tJCI/2, This is an 
interesting result for it implies that the PBCS state (j^) 
can be approximated to a great accuracy with a number 
of particle and hole states of order y/Q/2. We expect 
this result to hold for generic PBCS states. 

The gaussian decay of the weights we of the density 
matrix offers an ideal situation for the application of the 
DMRG. The DMRG works very well in the cases where 
the weights decay exponentially ( see jl3| for other types 
of decays). The p-h breaking allows for a smooth evo- 
lution of the system from a few electron regime into a 
superconducting one. 

Before we introduce the DMRG it is convenient to per- 
form the following canonical transformation, 
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j = i,...,n/2 
j = n/2 + i, ft 
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where the operators at a ( resp. &t a ) create particles ( 
resp. holes) acting on the Fermi sea (||). Choosing the 
chemical potential jj, as 



M=-(ft + l-A) 

the Hamiltonian (Q) has the p-h symmetry a,j 
and it can be written as, 



(6) 



H/d = K A + K B - \{A^A + B^B + AB + A^B^) (7) 
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where K B and B can be obtained from K A and A by 
the p-h transformation a^ a <-> bi^ a . We have substracted 
in (^) the constant term — (^) 2 , which is the energy 
of Fermi sea \FS), so that the lowest eigenvalue of H 
gives directly the condensation energy for even grains 
E§ = (ip\H\ip) - (FS\H\FS). For odd grains the level 
located at the Fermi sea is blocked and the condensation 
energy can be computed following the same steps as in 
the even case. From now on we shall concentrate on the 
latter case. We are now ready to apply the DMRG to 
find the GS of the Hamiltonian (]?]). At half filling this 
state takes the generic form 



\i>) = 



(10) 
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where the particle state \a) p must have a number of par- 
ticles N p equal to the number of holes Nh of \P)h for 
the GS state @ to be non vanishing. The DMRG is 
an algorithm that gives an optimal choice for the set of 
particle and hole states entering in (|l0|). This set is con- 
structed in successive steps starting from small grains. 
We begin with a system with n = 4 energy levels, which 
are chosen as the closest two particle and hole states near 
€f- This system can be represented as • • oo, where • 
stands for a particle level, while o stands for a hole one. 
The Fermi energy lies in between the »'s and the o's. The 
next step is to look for the GS of the Hamiltonian (Q) for 
n = 4 in the sector N p = Nh- From the knowledge of 
ipa. p w e define the reduce density matrix for the particle 
subspace 
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The p-h symmetry implies that the corresponding den- 
sity matrix in the hole subspace coincides with (|Tl|). The 
particle states which contribute the most to the GS ( |To| ) 
are the eigenvectors of the density matrix p A with high- 
est eigenvalues w p 0. For the system • • oo we can work 
with all the eigenstates of p A , but in general we shall 
only be able to keep the m most probable ones. With 
the information gained previously one builds the system 
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with SI = 6. The general rule is to build the system with 
D, = 2(n + 1) levels out of the system with SI = In. This 
is achieved by constructing the system with fl = 2(n+l) 
as a superblock of the form uA n B n o, where A n (resp. 
B n ) is the block which gives an effective description of 
the lowest n particle ( resp. hole) levels in terms of m 
states, while • and o represent the (n + l) th particle and 
hole levels added to enlarge the system size. 

The Hamiltonian H,abo of the superblock •A n B n o is 
given by 

H»abo — Ha + Hb + H, + H 
+ Hab + H,a + Hao + H,b + Hbo + H, Q (12) 



H A 
H. 

Hab 
H,a 

Hao 

H, 



XA\A n 



-iai 



l,er a "+l. cr 

-X(A n B n + h.c.) 
-X(A n al l+1 + h.c.) 
-X(A n b n+ i + h.c.) 
-\(a n+ ib n+ i + h.c.) 



Xa n+1 a n+ i 



(13) 
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and A n , K A are the operators A and 
and (0) but with SI/ 2 replaced by 



where a\ 

K A defined in eqs. 

n. K^ B \bi and B n have similar definitions. The terms 
Hb, H , H,b and Hbo can be derived from ( |l3| ) by the p- 
h transformation A n <-> B n , a, «-> h. The splitting ( |l2| ) of 
the superblock Hamiltonian H,ab<> recalls the one used 
in the momentum space DMRG 13] . However the latter 
reference uses a finite system algorithm which does not 
exploits the p-h symmetry. The DMRG provides a many 
body description of the blocks A n and B n , which means 
that the operators acting on these blocks are represented 
bymxm matrices. In our case the operators that we need 
to keep track are [A n ], [A^An] and [o'flj,,]. Given these 
operators we can construct the superblock Hamiltonian 
( |l3[ ) and look for the GS in the sector N p = Nh using the 
Lanczos method. The dimension of this Hubert space ( 
dim7iQ im ) is smaller than 4m 2 , for the constraint N p — 
Nh eliminates the states away from half filling. dim7iQ im 
is usually much smaller than the exact dimension of the 
Hilbcrt space of states with SI levels at half filling which 
is given by the combinatorial number Cq ,n/2- 

Given the GS of the superblock we obtain, using 
eq.(^), the density matrix of the particle system »A n 
and diagonalize it keeping the m most probable states 
with weight w p . The error of the truncation is measured 
by 1 — P m (P m = Yl™=i w p)- The latter states form a new 
basis of 'An denoted as A n +i and they give an effective 
description of the particle subspace with n + 1 levels. 
The hole block B n +\ is a mirror image of the particle 
block An+\. The DMRG proposed above is an infinite 
system algorithm, which is sufficient to study moderate 
system sizes (TV < 400). A way to improve the numerical 
accuracy of the infinite system method is to choose an 
effective value of the coupling constant A„(bulk) at the 



n th DMRG step in such a way that the value of the bulk 
gap is the one of the final system with coupling constant 
A = A (bare). This is guaranteed by the equation 
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A„(bulk) 



2(n+ 1) 
SI 



sinh 



1 



A (bare) 



(14) 



where SI is the final number of levels to reach and 2(n+ 1) 
is the number of levels at each step. Let us now present 
our results. A DMRG calculation for SI — 24 and m = 60 
agrees with the exact Lanczos condensation energy ob- 
tained in in the first 9 digits. The largest DMRG su- 
perblock matrix involved in the calculation is 3066 to be 
compared with the Lanczos matrix of dimension 2704156. 
For SI < 400 and m = 60 the condensation energy is com- 
puted with a relative error less than 10~ 4 . 

The numerical improvement achieved by the use of the 
effective coupling constant A„(bulk) defined in ( |l4[ ) as 
compared with the use of A = A(bare) in the DMRG 
steps is illustrated in table 1. 



m 


££(bare)/d 


£^'(bulk)/d 


1 - P m 


dimWioo.m 


50 


-40.44623 


-40.50014 


2.0 x 10" 


-y 


2108 


70 


-40.48878 


-40.50068 


7.1 x 10" 


n 


3622 


90 


-40.49815 


-40.50074 


1.1 x 10- 


n 


6306 


110 


-40.49983 


-40.50075 


1.5 x 10- 


12 


9720 



Table 1. GS condensation energy Eq for SI = 100 and 
A = 0.4 computed using the effective coupling constant 
and the bare one. 1 — P m is the truncation error of the 
last iteration and dim7Yioo,m is the largest dimension of 
the superblock. 

Let us next consider the crossover between the weak 
and the strong coupling regimes. The PBCS results of 
ref. JllJ suggest a sharp crossover between these two dif- 
ferent regimes at characteristic level spacings dp = 0.5A 
and di = 0.25A. For d < the condensation energy is 
an extensive quantity (~ 1/d) corresponding to a BCS- 
like behavior, while for d > df this energy is an inten- 
sive quantity (almost independent of d) fllJl . In figure 1 
we plot the DMRG results together with those of refer- 
ence |ll}| . From this comparison it is apparent that the 
DMRG gives significant lower energies than the PBCS 
ansatz. The DMRG is a variational method and in the 
region under study we expect our results to coincide with 
the exact ones with a relative error less than lO -4 . Fig. 
1 also shows that the crossover between the BCS and the 
fluctuation dominated (f.d.) regimes takes place in a re- 
gion which is wider than the one predicted by the PBCS 
approach and that there is no signs of critical level spac- 
ings df . A more quantitative characterization of this 
crossover is obtained by fitting the DMRG results to the 
formula ( see Fig. 1) 

^ b c /A = - Cl /ln(l + C2 f) + C3 + C 4f (15) 
ci = 1.48, c 2 = 3.05, c 3 = -1.98, c 4 = 0.08, {b = 0) 
ci = 0.36, c 2 = 0.86, c 3 = -1.95, c 4 = 0.16, {b = 1) 
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FIG. 1. GS condensation energies E%(b = 0, 1) as a func- 
tion of d/A for A = 0.224.. £1 ranges from 22 ( resp. 23) up 
to 400 ( resp. 401) for even ( resp. odd) grains and m — 60. 
The PBCS results are those of ref. 0. 

which interpolates between the bulk-BCS like behavior, 
given by E§ = -f^A 2 /d ( d/A « 1), and the f.d. 
regime (d/A >> 1) characterized by logarithmic correc- 
tions ||. Indeed from ( |i~5| ) we get C2/C1 = 2.06 for 6 = 
and C2/C1 = 2.4 for 6=1 which are both close to the 
bulk value given by 2 ||. 

The previous results are consistent with the probabil- 
ities of the to states kept by the DMRG as a function 
of the number of particles (resp. holes) of the states 
\a) p ( resp. \ft)h) appearing in (H), which we denote 
as Np^h- Since the DMRG keeps several states with the 
same value of particles and holes, it is appropriate to sum 
the probabilities of all the states with the same value of 
Np—h, which we shall denote as w(N p -h)- Each value of 
w(Np-h) is dominated by a single state which contributes 
the most. In Fig. 2 we plot w(N p _h) in the case A = 0.224 
and several values of O, which correspond to those plot- 
ted in Fig.l. The rapid decay of the weights recalls the 
gaussian decay of the eigenvalues of the density matrix 
of the PBCS (§). 

For Q = 22, 100, 180, 270 and 400 the most probable 
states have N p -h = 0,0,1,1,2 respectively, while the 
next two most probable states have occupation numbers 
liVp-zji 1|. As f2 increases the most probable state moves 
to higher values of N p -h, becoming eventually commen- 
surable with Q in the extreme BCS regime. It may seem 
from these results that the p-h DMRG is not capable to 
describe the bulk-BCS regime. This is not so for we have 
indeed observed this regime for higher values of A. 

In summary, we have proposed in this letter a new ver- 
sion of the DMRG which is suitable to study the GS prop- 
erties of the BCS pairing Hamiltonian of the ultrasmall 
superconducting grains. We believe that the p-h DMRG 
formulated in this work can be applied to a more gen- 
eral variety of fermionic systems in Condensed Matter, 
Atomic, Molecular |l4| and Nuclear Physics. Performing 
a Hartree-Fock (HF) transformation we can look for the 
best single particle basis to begin the DMRG procedure. 




FIG. 2. DMRG weights for m = 60 and A = 0.224 as a 
function of the number of particle states for different values 
of SI. 

The main limitation concerns again the amount of im- 
portant states which in general will grow with the size of 
the system to some power. For example in the pairing 
problem it grows with the square root of the size which 
allows us to study large systems. 

The comparison of the DMRG and the PBCS results 
p"T[ shows no signs of critical level spacings separat- 
ing qualitative different regimes. We rather observe a 
smooth logarithmic-like crossover which contradicts the 
sharp crossover predicted by the PBCS ansatz. This lat- 
ter feature is an artifact of the PBCS method which is 
unable to capture the true nature of the crossover. In 
summary the fluctuations seem to play a major role for 
no too small grains. 
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